# Figure 3

library(openxlsx)
library(cowplot)
library(ggpubr)

setwd('/Users/drschaef/Google Drive/Projects/Adriana/Papers/Homophily/Output')
setwd('/Users/drschaef/My Drive/Projects/Adriana/Papers/Homophily/Output')

rawCoeff <- read.xlsx('/Users/drschaef/Google Drive/Projects/Adriana/Papers/Homophily/Output/siena output 231128.xlsx', sheet='toPlot1')
rawCoeff <- sapply( rawCoeff, as.numeric )

colors3 <- c('cornflowerblue','darkblue','red')
shapes3 <- c(16,17,15)

# ERGM coefficient plot
plotCoeffs3 <- function(plotDat,xRange) {
  nEffects <- dim(plotDat)[1]
  CIwidth <- unlist(plotDat[,c(2,4,6)]) * 1.96
  CIlow <- unlist(plotDat[,c(1,3,5)]) - CIwidth
  CIhigh <- unlist(plotDat[,c(1,3,5)]) + CIwidth
  #par(mar=c(3,15,3,1))
  plot(unlist(plotDat[,c(1,3,5)]), c((nEffects:1)+.2, (nEffects:1)+0, (nEffects:1)-.2), type='p', pch=rep(shapes3, each=nEffects), 
       cex=1.3, col=rep(colors3, each=nEffects), yaxt='n', ylab='', xlab='', xlim=xRange, cex.lab=cexSize, cex.axis=cexSize, ylim=c(.5,2.5))
  abline(v=0, col='gray50', lty=2)
  segments(y0=c((nEffects:1)+.2, (nEffects:1)+0, (nEffects:1)-.2), y1=c((nEffects:1)+.2, (nEffects:1)+0, (nEffects:1)-.2), x0=CIlow, x1=CIhigh, col=rep(colors3, each=nEffects), lwd=1.75)
  title(xlab = "Estimate", line = 2.5, cex.lab=cexSize)            # Add x-axis text
}

modelNames <- c('Peer Influence \non Change','Peer Influence \non Stability')

# side-by-side version of Coefficient plot
xRange1 <- c(-.5,1.3)
xRange2 <- c(-.5,1.3)
cexSize <- .75
#tiff('fig4_coeff.tif', width=12, height=5, units='cm', res=300)
pdf('fig_coeff_231128.pdf', width=4.5, height=2)
hSetting <- length(modelNames)
layout(matrix(c(1:3), ncol=3), widths=c(1,2,2))
par(mar=c(4,0,2,.5))
plot(x=rep(1, hSetting),y=1:hSetting,xaxt='n',yaxt='n',xlab='',ylab='',pch='',bty='n',xlim=c(0,2), bg='white', ylim=c(.5,2.5))
text(rep(1, hSetting), (1:hSetting)+.2, rev(modelNames), cex = cexSize, pos=1, adj=c(0,0), ylim=c(.5,2.5))
#par(mar=c(4,0,3,1))
plotCoeffs3(rawCoeff[2:3,c(2,3,5,6,8,9)], xRange1)
mtext('MW',side = 3, at = sum(xRange1)/2, line = .5, cex=cexSize*.75) 
#legend("topleft", ncol=1, legend=c('Year 1', 'Year 2'), col=rev(colors2), pch=c(16,16), pt.bg= colors2, pt.cex=1, inset=.01, cex=cexSize)
plotCoeffs3(rawCoeff[4:5,c(2,3,5,6,8,9)], xRange2)
mtext('SW',side = 3, at = sum(xRange2)/2, line = .5, cex=cexSize*.75) 
dev.off()